! This module is an add on to ElastoDyn to allow output of blade structural data at each blade node when BeamDyn is not used
!
!  Copyright 2016   Envision Energy
!

MODULE ElastoDyn_AllBldNdOuts_IO

   USE NWTC_Library
   USE ElastoDyn_Types

   IMPLICIT NONE

   PRIVATE


   PUBLIC   :: AllBldNdOuts_InitOut
   PUBLIC   :: Calc_WriteAllBldNdOutput
   PUBLIC   :: AllBldNdOuts_SetParameters


      ! Parameters related to output length (number of characters allowed in the output data headers):

!   INTEGER(IntKi), PARAMETER      :: OutStrLenM1 = ChanLen-6    ! The NREL allowed channel name length is usually 20.  We are making these of the form B#N###namesuffix


! ===================================================================================================
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary.
! ===================================================================================================
! This code was generated by Write_ChckOutLst.m at 14-Dec-2017 10:34:30.


     ! Indices for computing output channels:
     ! NOTES:
     !    (1) These parameters are in the order stored in "OutListParameters.xlsx"
     !    (2) Array AllOuts() must be dimensioned to the value of the largest output parameter


     !  Local Span Motions:

   INTEGER(IntKi), PARAMETER      :: BldNd_ALx   =  1
   INTEGER(IntKi), PARAMETER      :: BldNd_ALy   =  2
   INTEGER(IntKi), PARAMETER      :: BldNd_ALz   =  3
   INTEGER(IntKi), PARAMETER      :: BldNd_TDx   =  4
   INTEGER(IntKi), PARAMETER      :: BldNd_TDy   =  5
   INTEGER(IntKi), PARAMETER      :: BldNd_TDz   =  6
   INTEGER(IntKi), PARAMETER      :: BldNd_RDx   =  7
   INTEGER(IntKi), PARAMETER      :: BldNd_RDy   =  8
   INTEGER(IntKi), PARAMETER      :: BldNd_RDz   =  9


     !  Local Span Loads:

   INTEGER(IntKi), PARAMETER      :: BldNd_MLx   = 10
   INTEGER(IntKi), PARAMETER      :: BldNd_MLy   = 11
   INTEGER(IntKi), PARAMETER      :: BldNd_MLz   = 12
   INTEGER(IntKi), PARAMETER      :: BldNd_FLx   = 13
   INTEGER(IntKi), PARAMETER      :: BldNd_FLy   = 14
   INTEGER(IntKi), PARAMETER      :: BldNd_FLz   = 15
   INTEGER(IntKi), PARAMETER      :: BldNd_MLxNT = 16
   INTEGER(IntKi), PARAMETER      :: BldNd_MlyNT = 17
   INTEGER(IntKi), PARAMETER      :: BldNd_FLxNT = 18
   INTEGER(IntKi), PARAMETER      :: BldNd_FlyNT = 19


     ! The maximum number of output channels which can be output by the code.
   INTEGER(IntKi), PARAMETER, PUBLIC      :: BldNd_MaxOutPts   = 19

!End of code generated by Matlab script
! ===================================================================================================


CONTAINS
!----------------------------------------------------------------------------------------------------------------------------------

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine populates the headers with the blade node outputs.  The iteration cycle is blade:node:channel (channel iterated
!! fastest).  If this iteration order is changed, it should be changed in the Calc_WriteAllBldNdOutput routine as well.
SUBROUTINE AllBldNdOuts_InitOut( InitOut, p, ErrStat, ErrMsg )


   TYPE(ED_InitOutputType),      INTENT(INOUT)  :: InitOut                          ! output data
   TYPE(ED_ParameterType),       INTENT(IN   )  :: p                                ! The module parameters
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat                          ! The error status code
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg                           ! The error message, if an error occurred

   INTEGER(IntKi)                               :: INDX                             ! Index count within WriteOutput
   INTEGER(IntKi)                               :: IdxBlade                         ! Counter to which blade we are on
   INTEGER(IntKi)                               :: IdxNode                          ! Counter to the blade node we ae on
   INTEGER(IntKi)                               :: IdxChan                          ! Counter to the channel we are outputting.
   CHARACTER(62)                                :: ChanPrefix                       ! Name prefix (B#N###)
   CHARACTER(6)                                 :: TmpChar                          ! Temporary char array to hold the node digits (3 places only!!!!)
   CHARACTER(*), PARAMETER                      :: RoutineName = ('AllBldNdOuts_InitOut')


         ! Initialize some things
      ErrMsg = ''
      ErrStat = ErrID_None

         ! First set a counter so we know where in the output array we are in
         ! NOTE: we populate invalid names as well (some names are not valid outputs for certain configurations).  That means we will have zeros in those values.
      INDX = p%NumOuts + 1       ! p%NumOuts is the number of outputs from the normal ElastoDyn output.  The WriteOutput array is sized to p%NumOuts + num(AllBldNdOuts)

      ! Populate the header and unit lines for all blades and nodes
#ifndef OUTPUT_CHANNEL_NAMES_AS_DISTANCE
         !   ! Warn if we will run into issues with more than 999 nodes.
         IF (p%BldNodes > 999 ) CALL SetErrStat(ErrID_Severe,'More than 999 blade nodes in use.  Output channel headers will not '// &
               'correctly reflect blade stations beyond 999. Modifications to the variable ChanLen in FAST are required.',ErrStat,ErrMsg,RoutineName)

            ! Populate the header an unit lines for all blades and nodes
            ! Loop through all the outputs we requested here:
         DO IdxChan=1,p%BldNd_NumOuts

            DO IdxBlade=1,p%BldNd_BladesOut
               DO IdxNode=1,p%BldNodes

                     ! 3 digit node number
                  WRITE (TmpChar,'(I3.3)')  IdxNode
                  ChanPrefix = 'B' // TRIM(Num2LStr(IdxBlade)) // 'N' // TRIM(TmpChar)

                     ! Now write to the header
                  InitOut%WriteOutputHdr(INDX) = TRIM(ChanPrefix) // TRIM(p%BldNd_OutParam(IdxChan)%Name)
                  InitOut%WriteOutputUnt(INDX) = p%BldNd_OutParam(IdxChan)%Units

                     ! Increment the index to the Header arrays
                  INDX = INDX + 1
               ENDDO

            END DO
         END DO

#else
      ! output format the name of the channel with the distance in cm from the root of the blade instead of by node number.

            ! Loop through all the outputs we requested here:
         DO IdxChan=1,p%BldNd_NumOuts
            DO IdxBlade=1,p%BldNd_BladesOut
               DO IdxNode=1,p%BldNodes

                     ! Node defined by distance along blade
                  WRITE (TmpChar,'(I6.6)')  NINT( 1000.0_ReKi * p%RNodes( IdxNode ) )
                  ChanPrefix = 'B' // TRIM(Num2LStr(IdxBlade)) // '_Z' // TRIM(TmpChar) //'_'


                     ! Now write to the header
                  InitOut%WriteOutputHdr(INDX) = TRIM(ChanPrefix) // TRIM(p%BldNd_OutParam(IdxChan)%Name)
                  InitOut%WriteOutputUnt(INDX) = p%BldNd_OutParam(IdxChan)%Units

                     ! Increment the index to the Header arrays
                  INDX = INDX + 1
               ENDDO

            ENDDO
         ENDDO

#endif


END SUBROUTINE AllBldNdOuts_InitOut

!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine populates the headers with the blade node outputs.  The iteration cycle is blade:node:channel (channel iterated
!! fastest).  If this iteration order is changed, it should be changed in the Calc_WriteAllBldNdOutput routine as well.
SUBROUTINE Calc_WriteAllBldNdOutput( p, u, m, y, LinAccES, ErrStat, ErrMsg )
   TYPE(ED_ParameterType),    INTENT(IN   )  :: p                                 ! The module parameters
   TYPE(ED_InputType),        INTENT(IN   )  :: u                                 ! inputs
   TYPE(ED_MiscVarType),      INTENT(INOUT)  :: m                                 ! misc variables
   TYPE(ED_OutputType),       INTENT(INOUT)  :: y                                 ! outputs
   REAL(ReKi),                INTENT(IN   )  :: LinAccES(:,0:,:)                  ! Total linear acceleration of a point on a   blade (point S) in the inertia frame (body E for earth). NOTE: zero index start.
   INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat                           ! The error status code
   CHARACTER(*),              INTENT(  OUT)  :: ErrMsg                            ! The error message, if an error occurred

      ! local variables
 
   INTEGER(IntKi)                               :: OutIdx                           ! Index count within WriteOutput
   INTEGER(IntKi)                               :: IdxBlade                         ! Counter to which blade we are on
   INTEGER(IntKi)                               :: IdxNode                          ! Counter to the blade node we ae on
   INTEGER(IntKi)                               :: J                                ! Generic counter for moment and force summation
   INTEGER(IntKi)                               :: IdxChan                          ! Counter to the channel we are outputting.
   CHARACTER(5)                                 :: ChanPrefix                       ! Name prefix ( B#N### or B#D#### )
   CHARACTER(2)                                 :: TmpChar                          ! Temporary char array to hold the node digits (2 places only!!!!)
   CHARACTER(*), PARAMETER                      :: RoutineName = 'Calc_WriteAllBldNdOutput'
!   REAL(ReKi)                                   :: ct, st ! cosine, sine of theta
!   REAL(ReKi)                                   :: cp, sp ! cosine, sine of phi
!   REAL(ReKi)                                   :: Tmp(3)
   REAL(ReKi)                                   :: OutVal                           ! Temporary variable to hold the value to output to the channel.

      ! Variables used in the CalcOutput routine that are needed here for coordinate transforms
   REAL(R8Ki)                   :: rSPS         (3)                                 ! Position vector from the undeflected blade node (point S prime) to the deflected node (point S)
   REAL(R8Ki)                   :: TmpVec       (3)                                 ! A temporary vector used in various computations.
   REAL(R8Ki)                   :: TmpVec2      (3)                                 ! A temporary vector.

 
         ! Initialize some things
      ErrMsg = ''
      ErrStat = ErrID_None


         ! Populate the header an unit lines for all blades and nodes
         ! First set a counter so we know where in the output array we are in
      OutIdx = p%NumOuts + 1       ! p%NumOuts is the number of outputs from the normal ElastoDyn output.  The WriteOutput array is sized to p%NumOuts + num(AllBldNdOuts)

         ! Loop through all the outputs we requested here:
      DO IdxChan=1,p%BldNd_NumOuts
      
            SELECT CASE( p%BldNd_OutParam(IdxChan)%Indx )      ! Indx contains the information on what channel should be output
            CASE (0) ! Invalid channel
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = 0.0_ReKi
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_ALx )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( LinAccES(:,IdxNode,IdxBlade), m%CoordSys%n1(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_ALy )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( LinAccES(:,IdxNode,IdxBlade), m%CoordSys%n2(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_ALz )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( LinAccES(:,IdxNode,IdxBlade), m%CoordSys%n3(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_TDx )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     rSPS = m%RtHS%rS0S(:,IdxBlade,IdxNode) - p%RNodes(IdxNode)*m%CoordSys%j3(IdxBlade,:)
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( rSPS, m%CoordSys%j1(IdxBlade,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_TDy )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     rSPS = m%RtHS%rS0S(:,IdxBlade,IdxNode) - p%RNodes(IdxNode)*m%CoordSys%j3(IdxBlade,:)
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( rSPS, m%CoordSys%j2(IdxBlade,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_TDz )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     rSPS = m%RtHS%rS0S(:,IdxBlade,IdxNode) - p%RNodes(IdxNode)*m%CoordSys%j3(IdxBlade,:)
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( rSPS, m%CoordSys%j3(IdxBlade,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_RDx )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( m%RtHS%AngPosHM(:,IdxBlade,IdxNode), m%CoordSys%j1(IdxBlade,:) )*R2D
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_RDy )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( m%RtHS%AngPosHM(:,IdxBlade,IdxNode), m%CoordSys%j2(IdxBlade,:) )*R2D
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_RDz )      ! See note in ElastoDyn.f90
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     ! y%WriteOutput( OutIdx ) = DOT_PRODUCT( m%RtHS%AngPosHM(:,IdxBlade,IdxNode), m%CoordSys%j3(IdxBlade,:) )*R2D           ! this is always zero for FAST
                     y%WriteOutput( OutIdx ) = 0.0_Reki
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_FLx )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( FrcMGagB(), m%CoordSys%n1(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_FLy )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( FrcMGagB(), m%CoordSys%n2(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_FLz )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( FrcMGagB(), m%CoordSys%n3(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_MLx )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( MomMGagB(), m%CoordSys%n1(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_MLy )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( MomMGagB(), m%CoordSys%n2(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_MLz )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( MomMGagB(), m%CoordSys%n3(IdxBlade,IdxNode,:) )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

               ! Output blade loads in the blade coordinate system.
            CASE ( BldNd_FLxNT )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( FrcMGagB(), CoordSysNT1() )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_FLyNT )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( FrcMGagB(), CoordSysNT2() )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_MLxNT )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( MomMGagB(), CoordSysNT1() )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE ( BldNd_MLyNT )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = DOT_PRODUCT( MomMGagB(), CoordSysNT2() )
                     OutIdx = OutIdx + 1
                  END DO
               END DO

            CASE DEFAULT
               CALL SetErrStat( ErrID_Severe, "Coding error.  Output channel not properly set.",ErrStat,ErrMsg,RoutineName )
               DO IdxBlade=1,p%BldNd_BladesOut
                  DO IdxNode=1,p%BldNodes
                     y%WriteOutput( OutIdx ) = 0.0_ReKi
                     OutIdx = OutIdx + 1
                  END DO
               END DO

         END SELECT
      END DO
      
contains
function CoordSysNT1()
   REAL(ReKi)                   :: CoordSysNT1  (3)                                 ! A temporary matrix for removing the structural twist from the local output forces and moments, dimension 1
      ! removed
   CoordSysNT1(:)  =  p%CThetaS(IdxBlade,IdxNode)*m%CoordSys%n1(IdxBlade,IdxNode,:) + p%SThetaS(IdxBlade,IdxNode)*m%CoordSys%n2(IdxBlade,IdxNode,:)

end function CoordSysNT1

function CoordSysNT2()
   REAL(ReKi)                   :: CoordSysNT2  (3)                                 ! A temporary matrix for removing the structural twist from the local output forces and moments, dimension 2

   CoordSysNT2(:)  = -p%SThetaS(IdxBlade,IdxNode)*m%CoordSys%n1(IdxBlade,IdxNode,:) + p%CThetaS(IdxBlade,IdxNode)*m%CoordSys%n2(IdxBlade,IdxNode,:)

end function CoordSysNT2


function FrcMGagB()
   REAL(R8Ki)                   :: FrcMGagB     (3)                                 ! Total force  at the blade element   (body M) / blade strain gage location            (point S) due to the blade above the strain gage.

      ! Initialize FrcMGagB using the tip brake effects:
      
   FrcMGagB = m%RtHS%FSTipDrag(:,IdxBlade) - p%TipMass(IdxBlade)*( p%Gravity*m%CoordSys%z2 + LinAccES(:,p%TipNode,IdxBlade) )
      
      ! Integrate to find FrcMGagB and MomMGagB using all of the nodes / elements above the current strain gage location:
   DO J = ( IdxNode + 1 ),p%BldNodes ! Loop through blade nodes / elements above strain gage node
      
      TmpVec2  = m%RtHS%FSAero(:,IdxBlade,J) - p%MassB(IdxBlade,J)*( p%Gravity*m%CoordSys%z2 + LinAccES(:,J,IdxBlade) )  ! Portion of FrcMGagB associated with element J
      FrcMGagB = FrcMGagB + TmpVec2*p%DRNodes(J)
      
   ENDDO ! J - Blade nodes / elements above strain gage node
      
      ! Add the effects of 1/2 the strain gage element:
      ! NOTE: for the radius in this calculation, assume that there is no
      !   shortening effect (due to blade bending) within the element.  Thus,
      !   the moment arm for the force is 1/4 of p%DRNodes() and the element
      !   length is 1/2 of p%DRNodes().
      
   TmpVec2  = m%RtHS%FSAero(:,IdxBlade,IdxNode) - p%MassB(IdxBlade,IdxNode)* ( p%Gravity*m%CoordSys%z2 + LinAccES(:,IdxNode,IdxBlade) ) ! Portion of FrcMGagB associated with 1/2 of the strain gage element
   FrcMGagB = FrcMGagB + TmpVec2 * 0.5 * p%DRNodes(IdxNode)                                                    ! Portion of FrcMGagB associated with 1/2 of the strain gage element
   FrcMGagB = 0.001*FrcMGagB           ! Convert the local force to kN

end function FrcMGagB

function MomMGagB()
   REAL(ReKi)                   :: MomMGagB     (3)                                 ! Total moment at the blade element   (body M) / blade strain gage location            (point S) due to the blade above the strain gage.

      ! Initialize MomMGagB using the tip brake effects:
      
   TmpVec2  = m%RtHS%FSTipDrag(:,IdxBlade) - p%TipMass(IdxBlade)*( p%Gravity*m%CoordSys%z2 + LinAccES(:,p%TipNode,IdxBlade) ) ! Portion of FrcMGagB
   MomMGagB = CROSS_PRODUCT( m%RtHS%rS0S(:,IdxBlade,p%TipNode) - m%RtHS%rS0S(:,IdxBlade,IdxNode), TmpVec2 )
      
      ! Integrate to find FrcMGagB and MomMGagB using all of the nodes / elements above the current strain gage location:
   DO J = ( IdxNode + 1 ),p%BldNodes ! Loop through blade nodes / elements above strain gage node
      
      TmpVec2  = m%RtHS%FSAero(:,IdxBlade,J) - p%MassB(IdxBlade,J)*( p%Gravity*m%CoordSys%z2 + LinAccES(:,J,IdxBlade) )  ! Portion of FrcMGagB associated with element J
      TmpVec = CROSS_PRODUCT( m%RtHS%rS0S(:,IdxBlade,J) - m%RtHS%rS0S(:,IdxBlade,IdxNode), TmpVec2 )                     ! Portion of MomMGagB associated with element J
      MomMGagB = MomMGagB + ( TmpVec + m%RtHS%MMAero(:,IdxBlade,J) )*p%DRNodes(J)
      
   ENDDO ! J - Blade nodes / elements above strain gage node
      
      ! Add the effects of 1/2 the strain gage element:
      ! NOTE: for the radius in this calculation, assume that there is no
      !   shortening effect (due to blade bending) within the element.  Thus,
      !   the moment arm for the force is 1/4 of p%DRNodes() and the element
      !   length is 1/2 of p%DRNodes().
      
   TmpVec2  = m%RtHS%FSAero(:,IdxBlade,IdxNode) - p%MassB(IdxBlade,IdxNode)* ( p%Gravity*m%CoordSys%z2 + LinAccES(:,IdxNode,IdxBlade) ) ! Portion of FrcMGagB associated with 1/2 of the strain gage element
   TmpVec = CROSS_PRODUCT( ( 0.25_R8Ki*p%DRNodes(IdxNode) )*m%CoordSys%j3(IdxBlade,:), TmpVec2 )                              ! Portion of MomMGagB associated with 1/2 of the strain gage element
      
   MomMGagB = MomMGagB + ( TmpVec + m%RtHS%MMAero(:,IdxBlade,IdxNode) )* ( 0.5 *p%DRNodes(IdxNode) )
   MomMGagB = 0.001*MomMGagB           ! Convert the local moment to kN-m
   
end function MomMGagB

END SUBROUTINE Calc_WriteAllBldNdOutput


!----------------------------------------------------------------------------------------------------------------------------------
!> This routine validates and sets the parameters for the nodal outputs.
SUBROUTINE AllBldNdOuts_SetParameters( p, InputFileData, ErrStat, ErrMsg )
!..................................................................................................................................


      ! Passed variables:

   TYPE(ED_InputFile),           INTENT(IN   )  :: InputFileData    !< Data stored in the module's input file
   TYPE(ED_ParameterType),       INTENT(INOUT)  :: p                !< Parameters
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat          !< Error status of the operation
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg           !< Error message if ErrStat /= ErrID_None

   
      ! local variables
   integer(IntKi)                           :: k                                 ! Blade number
   integer(IntKi)                           :: j                                 ! node number
   character(*), parameter                  :: RoutineName = 'AllBldNdOuts_ValidateInputData'
   
   ErrStat = ErrID_None
   ErrMsg  = ""

      ! Check if the requested blades exist
   IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) THEN
      p%BldNd_BladesOut = 0_IntKi
   ELSE IF ((InputFileData%BldNd_BladesOut > p%NumBl) ) THEN
      CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be less than "//TRIM(Num2LStr(p%NumBl))//".", ErrStat, ErrMsg, RoutineName)
      p%BldNd_BladesOut = p%NumBl ! NOTE: we are forgiving and plateau to numBlades      
   ELSE
      p%BldNd_BladesOut = InputFileData%BldNd_BladesOut
   ENDIF


      ! Check if the requested blade nodes are valid
   ! InputFileData%BldNd_BlOutNd



      ! Set the parameter to store number of requested Blade Node output sets
   IF ( p%BD4Blades .and. InputFileData%BldNd_NumOuts > 0 ) THEN
      p%BldNd_BladesOut = 0_IntKi
      p%BldNd_NumOuts = 0_IntKi
      CALL SetErrStat( ErrID_Warn,' AllBldNdOuts option not available in ElastoDyn when BeamDyn is used.  Turning off these outputs.',ErrStat,ErrMsg,"SetPrimaryParameters" )
   ELSE
      p%BldNd_NumOuts = InputFileData%BldNd_NumOuts
   ENDIF
   if (p%BldNd_BladesOut==0) p%BldNd_NumOuts = 0

      ! Set the total number of outputs ( requested channel groups * number requested nodes * number requested blades )
   p%BldNd_TotNumOuts = p%BldNodes*p%BldNd_BladesOut*p%BldNd_NumOuts !p%BldNd_NumOuts * size(p%BldNd_BlOutNd) * size(p%BldNd_BladesOut)

!      ! Check if the blade node array to output is valid: p%BldNd_BlOutNd 
!      ! TODO: this value is not read in by the input file reading yet, so setting to all blade nodes
!      !        -- check if list handed in is of nodes that exist (not sure this is ever checked)
!      !        -- copy values over
!
!      ! Temporary workaround here:
!   ALLOCATE ( p%BldNd_BlOutNd(1:p%BldNodes) , STAT=ErrStat2 )
!   IF ( ErrStat2 /= 0_IntKi )  THEN
!      CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the ElastoDyn OutParam array.", ErrStat, ErrMsg, RoutineName )
!      RETURN
!   ENDIF
!   DO I=1,p%BldNodes          ! put all nodes in the list
!      p%BldNd_BlOutNd(i) = i
!   ENDDO


!      ! Check if the requested blades are actually in use: 
!      ! TODO: this value is not read in by the input file reading yet, so setting to all blades
!      !        -- check if list handed in is of blades that exist (not sure this is ever checked)
!      !        -- copy values over
!   ALLOCATE ( p%BldNd_BladesOut(1:p%NumBl), STAT=ErrStat2 )
!   IF ( ErrStat2 /= 0_IntKi )  THEN
!      CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the ElastoDyn OutParam array.", ErrStat, ErrMsg, RoutineName )
!      RETURN
!   ENDIF
!   DO I=1,p%NumBl        ! put all blades in the list
!      p%BldNd_BladesOut(i) = i
!   ENDDO

   if (p%BldNd_TotNumOuts > 0) then
      call BldNdOuts_SetOutParam(InputFileData%BldNd_OutList, p, ErrStat, ErrMsg )
         if (ErrStat >= AbortErrLev) return
   end if


END SUBROUTINE AllBldNdOuts_SetParameters
 

!**********************************************************************************************************************************
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these 
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary. 
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine checks to see if any requested output channel names (stored in the OutList(:)) are invalid. It returns a 
!! warning if any of the channels are not available outputs from the module.
!!  It assigns the settings for OutParam(:) (i.e, the index, name, and units of the output channels, WriteOutput(:)).
!!  the sign is set to 0 if the channel is invalid.
!! It sets assumes the value p%NumOuts has been set before this routine has been called, and it sets the values of p%OutParam here.
!! 
!! This routine was generated by Write_ChckOutLst.m using the parameters listed in OutListParameters.xlsx at 14-Dec-2017 10:34:30.
SUBROUTINE BldNdOuts_SetOutParam(BldNd_OutList, p, ErrStat, ErrMsg )
!..................................................................................................................................

   IMPLICIT                        NONE

      ! Passed variables

   CHARACTER(ChanLen),        INTENT(IN)     :: BldNd_OutList(:)                  !< The list out user-requested outputs
   TYPE(ED_ParameterType),    INTENT(INOUT)  :: p                                 !< The module parameters
   INTEGER(IntKi),            INTENT(OUT)    :: ErrStat                           !< The error status code
   CHARACTER(*),              INTENT(OUT)    :: ErrMsg                            !< The error message, if an error occurred

      ! Local variables

   INTEGER                      :: ErrStat2                                        ! temporary (local) error status
   INTEGER                      :: I                                               ! Generic loop-counting index
   INTEGER                      :: J                                               ! Generic loop-counting index
   INTEGER                      :: INDX                                            ! Index for valid arrays

   LOGICAL                      :: InvalidOutput(1:BldNd_MaxOutPts)                      ! This array determines if the output channel is valid for this configuration
   CHARACTER(ChanLen)           :: OutListTmp                                      ! A string to temporarily hold OutList(I)
   CHARACTER(*), PARAMETER      :: RoutineName = "BldNdOuts_SetOutParam"

   CHARACTER(OutStrLenM1), PARAMETER  :: ValidParamAry(42) =  (/  &   ! This lists the names of the allowed parameters, which must be sorted alphabetically
                               "ALX  ","ALY  ","ALZ  ","AX   ","AY   ","AZ   ","FLX  ","FLXNT","FLY  ","FLYNT", &
                               "FLZ  ","FLZNT","FX   ","FXL  ","FY   ","FYL  ","FZ   ","FZL  ","MLX  ","MLXNT", &
                               "MLY  ","MLYNT","MLZ  ","MLZNT","MX   ","MXL  ","MY   ","MYL  ","MZ   ","MZL  ", &
                               "RDX  ","RDY  ","RDZ  ","RX   ","RY   ","RZ   ","TDX  ","TDY  ","TDZ  ","UXB  ", &
                               "UYB  ","UZB  "/)
   INTEGER(IntKi), PARAMETER :: ParamIndxAry(42) =  (/ &                            ! This lists the index into AllOuts(:) of the allowed parameters ValidParamAry(:)
                                  BldNd_ALx ,   BldNd_ALy ,   BldNd_ALz ,   BldNd_ALx ,   BldNd_ALy ,   BldNd_ALz ,   BldNd_FLx , BldNd_FLxNT ,   BldNd_FLy , BldNd_FlyNT , &
                                  BldNd_FLz ,   BldNd_FLz ,   BldNd_FLx , BldNd_FLxNT ,   BldNd_FLy , BldNd_FlyNT ,   BldNd_FLz ,   BldNd_FLz ,   BldNd_MLx , BldNd_MLxNT , &
                                  BldNd_MLy , BldNd_MlyNT ,   BldNd_MLz ,   BldNd_MLz ,   BldNd_MLx , BldNd_MLxNT ,   BldNd_MLy , BldNd_MlyNT ,   BldNd_MLz ,   BldNd_MLz , &
                                  BldNd_RDx ,   BldNd_RDy ,   BldNd_RDz ,   BldNd_RDx ,   BldNd_RDy ,   BldNd_RDz ,   BldNd_TDx ,   BldNd_TDy ,   BldNd_TDz ,   BldNd_TDx , &
                                  BldNd_TDy ,   BldNd_TDz /)
   CHARACTER(ChanLen), PARAMETER :: ParamUnitsAry(42) =  (/  &  ! This lists the units corresponding to the allowed parameters
                               "(m/s^2)","(m/s^2)","(m/s^2)","(m/s^2)","(m/s^2)","(m/s^2)","(kN)   ","(kN)   ","(kN)   ","(kN)   ", &
                               "(kN)   ","(kN)   ","(kN)   ","(kN)   ","(kN)   ","(kN)   ","(kN)   ","(kN)   ","(kN-m) ","(kN-m) ", &
                               "(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ","(kN-m) ", &
                               "(deg)  ","(deg)  ","(deg)  ","(deg)  ","(deg)  ","(deg)  ","(m)    ","(m)    ","(m)    ","(m)    ", &
                               "(m)    ","(m)    "/)


      ! Initialize values
   ErrStat = ErrID_None
   ErrMsg = ""
   InvalidOutput = .FALSE.


!   ..... Developer must add checking for invalid inputs here: .....

!   ................. End of validity checking .................


   !-------------------------------------------------------------------------------------------------
   ! Allocate and set index, name, and units for the output channels
   ! If a selected output channel is not available in this module, set error flag.
   !-------------------------------------------------------------------------------------------------

   ALLOCATE ( p%BldNd_OutParam(1:p%BldNd_NumOuts) , STAT=ErrStat2 )
   IF ( ErrStat2 /= 0_IntKi )  THEN
      CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the ElastoDyn BldNd_OutParam array.", ErrStat, ErrMsg, RoutineName )
      RETURN
   ENDIF


      ! Set index, name, and units for all of the output channels.
      ! If a selected output channel is not available by this module set ErrStat = ErrID_Warn.

   DO I = 1,p%BldNd_NumOuts

      p%BldNd_OutParam(I)%Name  = BldNd_OutList(I)
      OutListTmp          = BldNd_OutList(I)
      p%BldNd_OutParam(I)%SignM = 1   ! this won't be used

      CALL Conv2UC( OutListTmp )    ! Convert OutListTmp to upper case


      Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )

      IF ( Indx > 0 ) THEN ! we found the channel name
         IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN  ! but, it isn't valid for these settings
            p%BldNd_OutParam(I)%Indx  = 0                 ! pick any valid channel (I just picked "Time=0" here because it's universal)
            p%BldNd_OutParam(I)%Units = "INVALID"
            p%BldNd_OutParam(I)%SignM = 0
         ELSE
            p%BldNd_OutParam(I)%Indx  = ParamIndxAry(Indx)
            p%BldNd_OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output
         END IF
      ELSE ! this channel isn't valid
         p%BldNd_OutParam(I)%Indx  = 0                    ! pick any valid channel (I just picked "Time=0" here because it's universal)
         p%BldNd_OutParam(I)%Units = "INVALID"
         p%BldNd_OutParam(I)%SignM = 0                    ! multiply all results by zero

         CALL SetErrStat(ErrID_Fatal, TRIM(p%BldNd_OutParam(I)%Name)//" is not an available output channel.",ErrStat,ErrMsg,RoutineName)
      END IF

   END DO

   RETURN
END SUBROUTINE BldNdOuts_SetOutParam
!----------------------------------------------------------------------------------------------------------------------------------
!End of code generated by Matlab script
!**********************************************************************************************************************************

END MODULE ElastoDyn_AllBldNdOuts_IO
